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Section  1 


Introduction 


OPTRAN  is  a  raytrace  code  which  evaluates  the  optical  quality  of  aircraft  tranparencies 
subjected  to  operation  load  conditions.  This  volume  describes  the  theoretical  background 
on  which  the  code  is  based. 

The  raytrace  optical  code  is  interfaced  to  finite  element  thermal  and  stress  codes  to 
permit  the  effects  of  operational  loads  to  be  modeled.  Thermal,  displacement,  and  stress 
field  definition  data  computed  by  the  finite  element  codes  are  input  to  the  optics  code. 
This  information  is  required  to  compute  the  orthotropic  indices  cr  refraction  throughout 
the  material  volume  of  the  aircraft  transparency.  This  computation  is  performed  at  each 
step  along  the  propagation  path  of  each  ray. 

The  optics  code  tracks  rays  of  various  wavelengths  through  the  transparency.  The 
deformed  geometry  generated  by  the  stress  analysis  is  used  to  determine  angles  of  reflec¬ 
tion  and  refraction  at  transparency  layer  boundaries.  Birefringent  indices  of  refraction 
are  computed  as  a  function  of  material,  temperature,  and  stress  state  at  the  refracting 
surfaces  and  within  the  transparency  material. 

Post-processing  graphics  codes  display  the  angular  deviation,  transmittance,  and  po¬ 
larization  effects  over  specified  regions  of  the  transparency.  Plots  of  displacement  vectors 
and  deformed  grids  can  be  also  generated. 

The  PATRAN  [1]  finite  element  pre-  and  post-processing  software  provides  the  com¬ 
mon  interface  between  the  thermal  and  stress  finite  element  codes  and  the  the  optical 
analysis  code.  PATRAN  is  a  software  product  of  the  PATRAN  Division  of  PDA  En¬ 
gineering,  Costa  Mesa,  CA.  PATRAN  software  is  available  for  many  computer  systems 
and  offers  device  drivers  for  many  interactive  graphics  terminals.  PATRAN  provides 
excellent  tools  for  defining  the  model  geometries  and  generating  graphic  displays  of  the 
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models.  It  is  widely  used,  well  supported,  and  interfaced  to  many  of  the  more  popu¬ 
lar  finite  element  codes,  thus  offering  the  opportunity  of  using  other  analysis  tools  in 
conjunction  with  OPTRAN. 

Isoparametric  cubic  hyperpatches  defined  by  the  same  mathematical  formulation  as 
those  used  by  PATRAN  [lj  (Chapter  37)  are  used  in  OPTRAN.  They  are  also  used 
to  model  the  deformed  transparency  layer  solid  geometry  and  map  the  temperature 
and  stress  parameters  within  the  transparency  material  layers.  The  temperature  (T) 
and  the  six  orthogonal  stress  parameters  (  avv  otl  tzv  ryz  txz  )  are  required  at 
each  incremental  step  along  the  optical  ray  trace  paths  to  compute  continuously  varying 
orthotropic  indices  of  refraction.  Displacements  are  critical  in  determining  reflected  and 
refracted  optical  ray  path  directions  at  layer  boundaries. 

Temperature,  pressure,  and  density  fields  in  the*surrounding  air  stream  can  also  be 
mapped  using  a  series  of  isoparametric  hyperpatches.  These  parameters  determine  index 
of  refraction  in  the  atmosphere. 

The  mathematical  formulation  of  hyperpatches  is  presented  in  detail  in  Chapter  37 
of  the  PATRAN  Plus  User’s  Manual  [l].  The  hyperpatch  formulation  is  that  of  a  cubic 
solid  (64  node)  isoparametric  finite  element.  Coordinates  and  other  data  parameters  are 
mapped  within  the  hyperpatch  with  the  same  parametric  equations. 

For  raytracing,  the  entrance  and  exit  surfaces  of  each  part  must  be  identified.  Surfaces 
must  be  numbered  sequentially  from  the  outside  of  the  aircraft  to  the  eye.  The  eye 
position  and  a  set  of  pilot  reference  axes  must  also  be  defined.  Rays  are  specified  by 
direction  angles  with  respect  to  the  pilot  coordinate  system.  Rays  can  also  be  specified 
indirectly  by  defining  a  mesh  of  nodes  over  the  first  entrance  surface.  The  3D  coordinates 
of  these  nodes  are  used  to  generate  ray  directions. 

The  OPTRAN  code  accounts  for  three-dimensional  orthotropic  optical  effects.  An 
orthotropic  index  of  refraction  ellipsoid  is  computed  as  a  function  of  the  stress  and  tem¬ 
perature  values.  Orthotropic  effects  can  result  from  either  orthotropic  material  properties 
or  from  form  birefringance  caused  by  stress  optic  effects. 

Optical  material  properties  include  orthotropic  indices  of  refraction,  orthotropic  tem¬ 
perature  coefficients  of  the  indices  of  refraction,  and  a  six  by  six  matrix  of  stress  optic 
coefficients.  Orthotropic  optical  material  properties  input  to  OPTRAN  are  defined  with 
respect  to  material  axes. 

The  material  axes  are  defined  with  respect  to  reference  axes.  The  reference  axes 
can  be  either  the  global  coordinate  axes  or  axes  defined  by  the  derivatives  of  the  global 
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coordinates  with  respect  to  the  geometric  hyperpatch  parametric  variables.  This  latter 
reference  axis  option  allows  the  optic  axis  orientation  to  vary  with  the  curvature  of  the 
transparency.  By  default  the  material  axes  are  aligned  with  the  reference  axes  and  the 
global  coordinate  axes  are  the  reference  axes. 

The  following  method  is  used  to  find  the  ray  from  a  known  target  location  that 
intersects  the  eye  point.  First  we  aim  a  ray  from  the  target  to  the  eye  and  trace  an 
actual  ray  until  the  ray  intersects  the  eye  plane.  Then  we  trace  differential  rays  (close 
to  the  original  ray),  differing  first  in  azimuth  and  then  in  elevation.  The  intersection 
of  these  rays  at  the  eye  plane  (XY  plane  in  pilot  coordinates)  are  used  to  genrate  a 
first-order  matrix  that  can  be  solved  to  give  a  correction  to  the  original  ray,  that  is  a 
new  ray  that  now  intersects  the  eye  point. 

An  extensive  list  of  variables  are  generated  at  each  intersection  point.  These  include 
the  hyperpatch  ID  and  face  number,  material  ID,  hyperpatch  parameters  (u,  v,  w ), 
corresponding  3D  coordinates  ( x ,  y,  z ),  the  direction  of  the  surface  normal,  the  direction 
of  the  refracted  ray,  a  reference  polarization  direction,  polarization  and  transmittance 
arrays,  and  auxiliary  variables  needed  to  generate  differential  rays.  A  detailed  surface- 
by-surface  list  of  this  information  can  be  generated  on  the  output  listing. 

There  are  three  output  files  generated  by  OPTRAN.  The  first  is  the  output  listing, 
which  contains  a  copy  of  the  input  parameters,  detailed  ray  trace  information,  error 
messages,  and  summary  tables.  The  second  is  a  PATRAN  nodal  results  file,  which 
c<  tains  14  columns  of  summary  data.  The  PATRAN  nodal  results  file  can  be  used  in 
PATRAN  to  produce  a  variety  of  3D  plots.  The  third  file  is  an  optical  results  file,  which 
contains  outline  segments  that  help  define  the  field  of  view  and  raytrace  information  on 
a  uniform  grid  of  azimuth/elevation  variables.  The  optical  results  file  is  used  to  generate 
two-dimensional  grid  distortion  plots,  angular  deviation  fields,  and  polarization  ellipses. 

1.1  Optical  Raytrace 

Rays  are  represented  as  straight  lines  in  space.  Rays  are  refracted  at  the  point  they  inter¬ 
sect  an  optical  surface.  The  two  operations  involved  in  raytracing,  therefore,  are  finding 
the  intersection  of  a  ray  with  the  surface  and  refracting  the  ray  at  the  surface.  For  para¬ 
metric  surfaces,  ray  intersection  is  an  iterative  procedure,  requiring  a  two-dimensional 
nonlinear  optimization.  Calculating  the  direction  of  propagation  of  the  refracted  ray  is 
also  an  iterative  process,  since  the  index  of  refractive  varies  with  direction. 

Finding  the  first  intersection  on  an  entrance  surface  requires  a  search  over  available 
patches.  For  an  intersection  point  to  be  valid,  the  parametric  variables  corresponding 
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to  the  intersection  point  must  lie  within  the  bounds  of  0  <  (  u,  v .  w  )  <  1.  Once 
the  entrance  hyperpatch  face  has  been  identified,  the  ray  may  be  traced  through  the 
part  without  additional  searching  because  the  PATRAN  geometry  file  identifies  adjacent 
hyperpatches.  After  a  ray  exits  a  part,  however,  another  search  of  a  list  of  patches  may 
be  required  to  find  the  entrance  into  the  next  part. 

Stress  and  temperature  affect  light  rays  in  two  ways.  First,  the  surfaces  of  the  trans¬ 
parency  deform  as  a  result  of  stress  and  temperature  changes,  so  that  the  geometry  of 
the  transparency  changes.  Second,  the  index  of  refraction  of  the  transparency  layers 
depends  on  both  temperature  and  stress.  OPTRAN  uses  the  isoparametric  surface  ge¬ 
ometry  definitions  from  PATRAN  to  locate  boundaries  for  ray  refraction  and  reflection. 
The  results  of  finite-element  heat  and  stress  programs  define  the  volumetric  temperature 
and  strt  s  states  of  the  transparency,  from  which  the  index  ellipsoid  is  calculated. 

At  an  interface  between  two  dielectric  materials,  a  plane  of  incidence  is  defined  by  the 
normal  vector  to  the  surface  and  the  direction  vector  of  the  incident  ray.  The  light  ray  is 
split  into  reflected  and  refracted  rays,  propagating  in  the  plane  of  incidence.  Snell’s  law 
determines  the  direction  of  propagation.  The  polarization  of  the  incident  ray,  defined  by 
the  electric  field  vector,  is  decomposed  into  components  parallel  and  perpendicular  to  the 
plane  of  incidence.  The  Fresnel  equations  determine  the  reflectance  and  transmittance 
rt  each  polarization  component. 

Within  a  birefringent  material,  the  electric  field  vector  must  be  decomposed  into 
components  parallel  to  the  principal  axes  of  the  dielectric  material,  as  determined  by  the 
dielectric  tensor.  The  components  propagate  with  slightly  different  phase  shifts  causing 
one  polarization  state  to  be  retarded  in  phase  with  respect  to  the  other. 


1.2  Optical  Analysis 

Rays  are  traced  from  the  outer  world  to  the  eye  pomt,  as  shown  in  Figure  1.1.  Azimuth 
and  elevation  angles  are  used  to  specify  the  actual  directi  >n  k  of  a  target  point  in  the 
outer  world.  This  is  the  dotted  line  drawn  from  the  object  to  the  eye  in  Figure  1.1.  The 
direction  of  the  exit  ray  k'  shows  the  apparent  direction  of  the  target,  as  shown  by  the 
dotted  line  from  the  image  to  the  eye.  The  difference  between  the  apparent  direction 
and  the  actual  direction  is  the  angular  deviation.  The  angular  deviation  6k  is  calculated 
from  6k  =  k'  -  k.  Deviation  causes  objects  to  be  seen  at  other  than  their  true  direction 
from  the  observer. 

In  OPTRAN,  2D  direction  coordinates  of  a  vector  are  defined  with  respect  to  the 
pilot  reference  system.  Note  that  direction  vectors  point  toward  the  eye,  and  that  zero 
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F*gure  1.1:  Angular  deviation  of  a  noint  • 

n  of  a  point  in  space  by  aircraft  transparency 


F,gUre  12  Diagram  illustrating 


the  definition  of  angular  coordinates 


angle  corresponds  to  a  direction  parallel  to  the  pilot  z-axis.  For  this  case  the  point  of 
regard  (object  or  image  point)  lies  along  the  negative  z-axis,  directly  in  front  of  the  pilot. 

The  magnitude  of  direction  coordinates  is  the  angle  <f>  in  degrees  between  the  vector 
and  the  pilot  z-axis,  as  shown  in  Figure  1.2.  The  orientation  of  the  angular  coordinates 
is  obtained  by  projecting  the  direction  vector  onto  the  xy-plane  and  finding  the  angle  9 
with  respect  to  the  x-axis,  as  shown  in  Figure  1.2.  Then  azimuth  and  elevation  direction 
angles  A x  and  Av  are  found  from 


Aa  —  <f>  cos  0 
Av  =  <j>  sin  6 

The  magnitude  of  direction  angles  extends  from  0  to  180°.  A  direction  angle  of  180° 
refers  to  a  point  directly  behind  the  pilot.  For  the  180°  direction  the  orientation  9  is 
indeterminate.  The  horizontal  or  x-component  of  direction  is  called  azimuth,  and  will 
be  positive  for  points  of  regard  to  the  right  of  the  pilot  and  negative  to  the  left.  The 
vertical  or  y-component  of  direction  is  called  elevation,  and  will  be  positive  for  points 
above  the  horizon  (xz-plane)  and  negative  below. 

Angular  deviation  is  the  most  direct  measure  of  optical  quality  of  a  transparency. 
One  method  of  demonstrating  angular  deviation  is  to  show  how  a  rectilinear  grid  in 
angle  coordinates  is  distorted  by  the  transparency.  The  data  from  a  square  array  of 
rays,  equally  spaced  in  azimuth  and  elevation,  can  be  displayed  in  single  cross-sections 
of  elevation  or  azimuth  error  or  as  deformed  grids,  emulating  the  typical  grid  board 
photograph. 

Grid  distortion  is  shown  in  Figure  1.3  (a).  The  undistorted  grid  is  shown  in  dashed 
lines  and  the  distorted  grid  by  solid  lines.  The  arrow  is  the  angular  deviation  for  a  single 
point  on  the  grid.  If  arrows  are  drawn  at  every  grid  point,  then  the  representation  shown 
in  Figure  1.3  (b)  is  obtained.  Finally,  by  removing  the  grids  we  can  produce  a  plot  of 
angular  deviation  vectors  as  shown  in  Figure  1.3  (c).  Both  grid  distortion  and  angular 
deviation  plots  can  be  produced  from  OPTRAN  optical  results  files. 


1.3  Polarization 

Characterizing  the  state  of  polarization  is  discussed  in  Born  and  Wolf  (3j  pp.  30-32, 
Hecht  (2]  pp.  321-326,  and  in  Azzam  and  Bashara  [9].  The  different  possible  states  of 
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Figure  1.3:  Using  a  grid  to  show  distortion  and  angular  deviation 

7 


polarization  of  monochromatic  light  can  be  represented  by  a  set  of  four  real  quantities, 
called  the  Stokes  parameters,  denoted  by  the  vector 

Sf  =  (  Sj  s2  s3  )  , 

each  component  of  which  has  the  dimensions  of  irradiance.  The  first  term,  s0,  gives  the 
total  irradiance  of  the  light  wave  and  therefore  is  always  positive.  The  next  three  terms 
give  the  difference  between  the  irradiance  of  three  different  sets  of  orthogonal  polarization 
states,  and  can  therefore  be  either  positive,  negative,  or  zero. 

Suppose  we  had  our  choice  of  ideal  linear  and  circular  polarizers.  Let  I0  be  the  total 
irradiance  of  the  wave  and  Asi  /_ 45,  A ,  It  be  the  intensities  transmitted  by  the 

corresponding  ideal  polarizer.  Then 


So  =  Iq  —  It  +  ly  —  As  +  I-  45  —  It  +  It 

Sj  —  /*  ~  ly 

Sz  =  As  —  1-45 

S3  —  Ir  ~  A 

where  gives  the  difference  between  the  irradiance  of  x  and  y  linear  polarization  states, 
s2  represents  the  preference  for  polarization  between  +45°  and  —45°,  and  s3  is  the  dif¬ 
ference  between  right-  and  left-circular  polarizations. 

For  unpolarized  light,  there  is  no  preference  for  any  particular  polarization  so  that 
Sj  =  s2  —  s3  =  0  and  the  Stokes  vector  has  the  simple  form 

sf  =  (  s0  0  0  0  )  . 

The  Stokes  parameters  of  a  totally  polarized  wave  satisfy  the  condition 

50  =  S1  +  S2  +  4 

The  general  case  can  be  treated  by  splitting  the  wave  into  two  components,  one  totally 
polarized  and  one  unpolarized.  Let 

sp  =  vsi +  4  +  4- 

Then 

=  (  (so  -  sP)  0  0  0  ) 

stp  -  {  SP  S 1  53  ) 
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Circular 


Elliptical 


Linear 


Figure  1.4:  Examples  of  polarization  states 


We  can  also  define  the  degree  of  polarization  Dp  as  the  ratio  of  the  irradiance  of  the 
totally  polarized  component  to  the  total  irradiance 


The  degree  of  polarization  varies  from  zero  for  unpolarized  light  to  unity  for  totally 
polarized  light,  with  intermediate  values  for  partially  polarized  light. 


Polarization  has  a  simple  geometric  interpretation  as  an  ellipse.  The  shape  and 
orientation  of  the  ellipse  relate  to  the  polarization  state,  as  shown  in  Figure  1.4.  The 
size  of  the  ellipse  relates  to  the  relative  brightness  of  the  polarized  component  of  the 
transmitted  ray.  The  shape  of  the  ellipse  may  be  characterized  by  two  angles,  the 
ellipticity  angle  c  and  the  orientation  angle  6 ,  as  shown  in  Figure  1.5.  The  relationship 
between  the  Stokes  parameters  and  the  polarization  ellipse  is 


Section  2 


Raytracing  Parametric  Surfaces 


Raytracing  is  the  process  of  following  a  ray  of  light  through  an  optical  system.  Rays  are 
represented  lines  in  space.  The  intersection  point  is  calculated  for  a  ray  with  a  surface. 
For  parametric  surfaces,  this  is  an  iterative  procedure,  requiring  a  two-dimensional  non¬ 
linear  optimization.  The  ray  is  next  refracted  or  reflected  at  the  surface,  and  a  new  ray 
direction  vector  is  calculated.  Then  the  raytrace  is  repeated  for  the  succeeding  surfaces 
until  the  image  surface  is  reached. 

This  chapter  presents  the  working  equations  necessary  to  find  the  intersection  between 
rays  and  parametric  surfaces. 

2.1  Overview 

Three  ray  transfer  mechanisms  have  been  devised  for  tracing  rays  through  optical  trans¬ 
parencies.  The  first  is  an  insertion  process.  Eligible  hyperpatches  for  ray  entry  are  those 
associated  with  the  entrance  surface.  Ray  intersections  are  tested  for  each  eligible  hy¬ 
perpatch  until  a  legal  intersection  is  found.  If  no  intersection  is  found,  the  action  taken 
depends  on  the  transmission  code  for  that  component.  If  the  transmission  code  is  zero, 
the  ray  is  assumed  to  be  blocked.  The  ray  trace  is  continued  to  the  next  component  if  the 
transmission  code  is  set.  The  second  mechanism  is  that  of  internal  propagation  in  which 
a  ray  is  traced  from  one  side  of  the  hyperpatch  to  another.  The  third  is  ray  extraction 
at  the  exit  face  of  a  hyperpatch.  The  adjoining  face  is  either  a  new  hyperpatch,  the  exit 
surface,  or  an  edge  face.  If  the  adjoining  face  is  a  new  hyperpatch,  the  local  parametric 
ray  coordinates  must  be  calculated  for  the  new  hyperpatch.  If  the  exit  surface  is  encoun¬ 
tered,  the  ray  is  refracted  into  air  and  traced  to  the  eye  plane  or  to  the  entrance  surface 
of  the  next  optical  component.  Finally,  if  an  edge  face  is  encountered,  the  ray  is  either 
reflected  from  the  interface  or  totally  absorbed  at  the  point  of  intersection. 
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2.2  Hyperpatches 


The  geometric  model  used  in  OPTRAN  is  the  hyperpatch  defined  by  PATRAN  (chapter 
37)  [lj.  A  hyperpatch  is  a  3-dimensional  mapping  of  the  parametric  unit  cube  in  (u,  v, 
w)  to  world  coordinates  (x,  y,  z). 

The  following  method  evaluates  a  tricubic  hyperpatch  at  (u,v,w)  to  obtain  (X,Y,Z) 
coordinates. 

Formal  Parameters: 


•  C  (Input)  =  Hyperpatch  coefficients.  For  each  coordinate  (X,Y,Z),  these  are  ar¬ 
ranged  in  a  4x4x4  matrix  of  the  form: 


c 

(0,0,0) 

c 

(0,1,0) 

cu 

c„ 

c 

(1,0,0) 

c 

(1,1,0) 

c„ 

10001 

c„ 

1800] 

c„ 

(0,0,0) 

ESE1 

Cuu 

cuv 

Cu 

(1,0,0) 

wmlm 

cuv 

■0061 

c„„ 

c 

(0,0,1) 

c 

(0,1,1) 

tarn 

KSSOl 

EC 

c 

(1,0,1) 

c 

(144) 

19 

worn 

EM 

mum 

Cuw 

IfiSs&l 

Cuv 

imau 

IQQQ1 

Cw 

KEfil 

Cuv 

■oxjxoi 

Cvw 

102291 

igkes 

Cvw 

■0001 

1800 

Cuw 

(0,0,0) 

Cuul 

(0,1,0) 

Cuvxo 

(0,0,0) 

CttVW 

(0,1,0) 

Cuw 

(1,0,0) 

Cuw 

(1,1,0) 

Cuvw 

(1,0,0) 

Cut,  ,0 

(1,1,0) 

cw 

wmsi 

cw 

(0,1,1) 

Cvw 

(0,0,1) 

Cvw 

iiaiu] 

cw 

■will 

Cw 

(1.1.1) 

Cvw 

(1,0,1) 

Cvw 

18001 

I0Q3P1 

CltO 

(0,1,1) 

Cuvw 

(0,0,1) 

ICES] 

Cuw 

(1,1,1) 

Cuvw 

(1,0,1) 

18001 

•  UVW  (Input)  =  Parametric  coordinates  on  the  interval  (0,1) 

•  XYZ  (Output)  =  Cartesian  coordinates  at  point  (u,v,w)  of  the  patch 


Method: 
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For  each  coordinate,  evaluate: 

Hu  =  MU 

H„  -  MV 

H,,  =  MW 

in  which: 

U*  =  (  u3  u2  u  1  ) 

V*  =  |  u3  u2  u  1  j 

Wf  =  (  w3  w2  w  1  ) 

C  is  the  4x4x4  matrix  of  patch  coefficients  for  coordinate  X,  and  M  is  a  4x4  matrix 
defining  the  coefficients  of  the  first  order  (cubic)  Hermite  polynomials: 


M  = 


2-301' 
-2  3  0  0 

1-210 
1-10  0 


=  (Hu)i  *  {Hv)j  *  {Hw)k 
l  —  i  +  4 [j  -  1)  +  16 [k  -  1) 


J  for  i=l-4,  j  =  l-4,  k= 1-4 


Note:  The  products  MU,  MV,  and  MW  are  the  same  for  all  three  coordinates.  The 
matrix  M  above  is  the  transpose  of  the  matrix  M  defined  in  the  PATRAN  theory  chapter. 


X  =  £ff,C,, 


i=i 


64 

y  =  64, 


i=i 


64 

Z  =  H1C1+ us 

»=i 


2.3  Parametric  Surfaces 

Parametric  surfaces  are  found  by  setting  one  of  the  parametric  variables  to  one  of  its 
limiting  values  (0  or  1).  This  yields  one  of  the  six  faces- of  the  hyperpatch.  Although 
any  one  of  the  six  faces  may  be  used  in  OPTRAN,  we  show  the  case  for  w  held  constant. 
This  parametric  surface  is  defined  by 

S(u,  v)  =  (  X(u,  v)  T(u,v)  Z(u,v)  ) 
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Surface  Normal 


At  any  point  (u,  v)  on  the  parametric  surface,  we  can  construct  a  vector  N  perpendicular 
to  the  patch  by  computing  the  cross  product  of  the  tangent  vectors  Su  and  S”. 

N  =  S"  x  S“ 


where 


S“ 

S” 


dS(u,v) 

du 

dS  (u,  u) 
dv 


2.4  Surface  Intersection 


A  general  discussion  of  geometric  modeling  based  on  parametric  cubic  surfaces  is  found  in 
Mortenson  [4].  Ray  tracing  of  parametric  surfaces  has  been  applied  to  studies  of  computer 
graphics  [5,6,7]. 

Given  the  parametric  surface 

S(u,u)  =  (  AT(u,v)  F(u,u)  Z(u,v)  ) 

and  a  ray  defined  by  the  point  p  and  the  unit  direction  k,  the  square  of  the  distance 
from  the  surface  to  the  ray  is 

<t>2  =  |f|2  =  |p'  -  p  -  <7k|’ 

where  F  is  the  defect  vector  and  p’  the  intersection  point, 

p'  =  S(u,u) 

and  the  distance  q  along  the  ray  to  the  point  of  intersection  is 

7  =  (P'  -  P)  ■  k. 

A  local  minimum  of  <t>2  =  0  corresponds  to  a  point  (u,  v)  where  the  ray  intersects  the 
surface.  All  components  of  the  defect  vector  F  must  vanish  for  an  intersection  point  to 
exist.  A  solution  can  be  obtained  by  adjusting  the  variables  u,  v,  and  q.  Those  points 
where  4> 3  is  a  minimum,  but  <f> 2  >  0,  indicate  that  the  ray  missed  the  surface  by  a  finite 
distance.  Thus  a  minimization  algorithm  will  still  converge  near  a  “silhouette  edge”  of 
the  surface,  but  will  give  a  non-zero  minimum. 
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2.5  Ray  Refraction 


Refraction  at  a  boundary  between  two  dielectric  media  is  determined  by  Snell’s  law,  that 

»(k  x  N)  =  n'(k'  x  N) 

where  N  is  the  surface  normal  vector  and  k  the  unit  vector  in  the  direction  of  propagation. 

We  can  obtain  another  representation  of  Snell’s  Law  by  examining  the  following 
vector  triple  product: 

n'(k'  x  N)  x  N  =  n(k  x  N)  x  N 
Using  a  standard  vector  identity,  we  can  write  this  as 

(N  •  N)n'k'  -  (N  ■  n'k')N  =  (N  •  N)nk  -  (N  •  nk)N 
Let 

g 2  =  N-N 
T  =  nk  •  N  =  ng  cos  6 
T'  =  n'k'  •  N  =  n'g  cos  O' 

These  definitions  allow  us  to  write 

n'k'  =  nk  +  7N 

where  ,  _  p 

1  =  - ’ 

r 

Now  let  us  address  the  problem  of  finding  P.  We  seek  a  solution  that  avoids  trigonometric 
functions. 

n'sintf'  =  nsinfl 
n'2(l  -  cos2  O')  =  n'2(l  -  cos2  0) 
g2n '2  -  T2  =  gW  -  T2 
r'2  =  r2  +  <72(n'2  -  n2) 

T'  =  ±yjT2  +g2{n'2  -nJ) 

If  the  argument  of  the  square  root  is  negative,  we  have  total  internal  reflection. 

The  choice  of  sign  can  be  made  as  follows: 


Refraction:  T'  and  T  have  the  same  sign 

Reflection:  P  and  f  change  signs 


For  reflection,  P  =  — P 
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2.6  Refraction  and  Reflection  Coefficients 


The  refraction  and  reflection  coefficients  determine  the  energy  which  is  transmitted  or 
reflected  at  a  dielectric  interface. 

The  calculation  of  these  coefficients  is  best  performed  in  a  local  coordinate  system 
defined  as  follows.  Let  z  be  the  normal  to  the  surface,  with  orthogonal  vectors  x  and  y 
lying  in  the  plane  of  the  surface,  and  let  x  be  perpendicular  to  the  plane  of  incidence, 
then 

x  =  k  x  z 
y  =  *  x  x 


If  the  propagation  vector  k0  is  parallel  to  the  surface  normal  (normal  incidence),  no  plane 
of  incidence  is  defined.  The  choice  of  x  is  then  arbitrary,  and  we  let  x  be  in  the  direction 
of  the  incident  polarization  reference  axis  d0. 


The  Maxwell  equations  for  reflection  and  refraction  at  the  boundary  between  two 
isotropic  materials  can  be  solved  for  ampltiude  reflection/refraction  coefficients, 

The  following  amplitude  reflection/ transmission  coefficients  may  be  obtained 


ti  = 


where  t  is  one  of  four  possible  combinations,  (reflecting  or  refracting)  and  (p-polarization 
or  s-polarization),  o  represents  the  entrance  medium  and  m  represents  the  exit  medium. 

Energy  reflection/transmission  coefficients  are  obtained  from 


T,  = 


I  umkm  •  z 
n0 k0  •  z 


For  reflection,  the  index  of  refraction  of  exit  and  entrance  media  are  the  same. 


S-Polarization 

For  an  incident  ray  polarized  perpendicular  to  the  plane  of  incidence  (s-polarization), 
the  electric  field  vector  e0  =  x  and  the  reflection  and  transmission  amplitude  coefficients 
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are 


2 n„ cos  9 

t,  =  - 2 - 

n0  cos  9  +  nm  cos  0' 

n0  cos  8  -  nm  cos  O' 

n0  cos  0  +  nm  cos  0' 

where  9  is  the  angle  of  incidence  and  01  the  angle  of  refraction.  For  this  case,  the 
amplitude  coefficients  simplify  to 

2r 

ta  ~  r  +  r 
r  -  r 

r*  -  r  +  r, 


P-Polari*ation 

For  an  incident  ray  polarized  parallel  to  the  plane  of  incidence  (p-polarization),  e0  =  y 
and 

2n0  cos  8 

t  —  - - - 

p  nm  cos  0  +  n0  cos  9' 

nm  cos  0  —  n0  cos  O' 

f  ~  ■ 

p  nm  cos  0  +  n„  cos  9' 


Normal  incidence 

For  normal  incidence,  there  is  no  plane  of  incidence  defined.  We  can  then  choose  x  =  e0. 
Then  the  equations  simplify  to 

Eh  _  2  n0 

E0  Tl0  -+■  Hm 

Ea  _  nm  —  n0 

E0  nm  +  n0 
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2.7  Mueller  Matrices 


In  a  general  polarizing  system,  the  output  polarization  state  is  related  to  the  input 
polarization  state  by 

5o  rnoo  m0i  rrio2  m03  /  5o 

Sj  _  rn,\ o  mu  mi2  m\z  Si 

s'2  rn2o  m2i  m22  m23  s2 

s3  /  m30  m31  m32  V  s3 

where  the  matrix  M  is  called  a  Mueller  matrix.  [9j 

Following  are  Mueller  matrices  for  common  optical  components.  The  Mueller  matrix 
for  a  polarizer  oriented  at  0°: 


110  0 
110  0 
0  0  0  0 
0  0  0  0 


Linear  retarder,  oriented  along  x-axis: 


1  0  0 

0  0  0 

0  0  cos  6 

0  —0  sin  6 


0  " 
0 

sin  6 
cos  6 


The  Mueller  matrix  for  reflection/ transmission  is  equivalent  to  a  partial  polarizer 
with  tranmission  factors  tx  tv 

tx  +  ty  tx  —  ty  0  0 

_  _  1  tX  ~  ty  tx  +  ty  0  0 

2  0  0  2y /tjTv  0 

000  2^/rj; . 

Rotation  matrix: 

1  0  0  0' 

0  cos  2a  sin  2a  0 

0  -  sin  2a  cos  2 a  0 
0  0  0  1. 

To  find  the  Mueller  matrix  of  a  device  oriented  at  an  azimuth  of  9: 

M  =  R(-0)MR(0) 
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2.8  Output  Data 


The  result  of  a  ray  trace  is  a  set  of  ray  data,  consisting  of  the  following  information  at 
intersection  points: 


NPAT 

geometric  hyperpatch  ID 

(u,  u,  w) 

parametric  coordinates 

p 

ray  position  vector  ( x ,  y,  z) 

n 

surface  normal  unit  vector 

k 

direction  of  wave  propagation 

d 

polarization  reference  vector 

N 

mean  index  of  refraction 

6N 

difference  in  refractive  index 

d 

distance  to  next  surface 

These  data  are  calculated  for  the  object  surface,  each  intermediate  surface  intersected 
by  the  ray,  and  finally  for  the  image  surface.  The  optical  path  to  the  next  surface  is  Nd. 
The  optical  retardance  is  6Nd. 
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Section  3 


Nonlinear  Least- Squares  Optimization 


The  surface  intersection  problem  defined  in  the  previous  section  is  solved  by  the  nonlinear 
least-squares  optimization  technique  presented  in  this  section. 

The  optimization  method  seeks  to  minimize  the  sum  of  squares  of  defect  terms,  and 
the  ideal  solution  is  obtained  if  each  individual  term  vanishes.  The  method  expands  the 
merit  function  in  a  Taylor  series  about  a  starting  point,  keeping  the  linear  and  quadratic 
terms  in  the  variables.  This  approximate  function  is  then  minimized,  and  a  solution 
vector  is  found.  This  vector  may  extend  outside  the  region  of  validity  for  the  Taylor 
series,  so  the  solution  vector  is  reduced  in  magnitude  or  damped  to  keep  changes  within 
the  region  of  validity. 

3.1  Mathematic  Preliminaries 

Merit  Function 

The  defect  functions  /,  are  functions  of  a  set  of  N  variables  (zj,  X2, . . . ,  X/y): 

fi  —  /i(zi>  x2i  •  ■  •  »  xn) 

fi  =  /j(lii^2i  ■  •  •  i^Af) 

/M  =  f\f(xl,x2,  ■  ■  ■  ixn) 

The  merit  function  is  of  the  type 

M 

<t>2  =  f' 

i  =  1 
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or 

<£2=ftf  =  f.f  =  ||f||* 

where  f  is  a  ( M  x  1)  vector  and  f*  is  the  (l  x  M )  transpose  of  f.  The  first  form  of  the 
expression  uses  the  notation  of  matrix  multiplication.  The  second  form  shows  a  vector 
dot  (or  inner)  product,  and  the  last  form  is  a  vector  norm  over  defect  space. 


Linear  Defect  Model 


Over  a  small  region  about  the  starting  point,  the  defects  may  be  approximated  by  a 
Taylor  series, 

f  =  ao  +  Ax. 

where  A  is  a  ( M  x  TV)  matrix  of  first  derivatives: 


and  x  are  changes  in  the  variables  from  the  starting  point. 


Gradient 

The  gradient  g  is  a  (N  x  1)  vector  given  by 


g  =  v<*2 


Its  components  are 


*  =  tel =  2  V'dZ  +  Aa^  +  "  +/MlJr  j 


then 

g  =  2A*f 

A  very  inefficient  method  for  finding  a  solution  point  is  to  search  for  a  minimum  along 
the  vector  -g.  This  is  known  as  the  gradient  search  method,  and  for  a  highly  nonlinear 
function  may  be  a  good  strategy.  Experience  has  shown,  however,  that  for  most  cases 
the  local  gradient  is  not  a  good  predictor  of  the  final  minimum. 
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3.2  Method  of  Least-Squares 

Using  the  linear  model  for  the  defects  allows  us  to  express  the  merit  function  as 
<t> 2  =  (a0  4-  Ax)  •  (a0  4  Ax)  =  a0  •  a0  +  2c0  •  x  +  x*Cx 


where 


c0  =  Afa0 
C  =  AfA 

Let  a;  represent  column  j  of  matrix  A.  The  matrix  C  is  a  symmetric  ( N  x  N)  matrix, 
whose  elements  can  be  written  as  a  sum  over  the  defects, 

M 

Cj  —  )  *  A,;(a,)o  —  a  j  •  ao 

•=i 

M 

Cjk  =  AjjiA,*  =  a  j  •  ak 

i=  1 

The  matrix  C  is  called  the  covariance  array.  The  gradient  is 

g  =  2Aff  =  2(c0  4-  Cx) 

The  minimum  of  </>2  is  obtained  by  setting  g  =  0  and  solving  for  x.  The  resulting  matrix 
equation 

c0  4-  Cx  =  0 

is  a  set  of  simultaneous  linear  equations  known  as  the  normal  equations  of  least-squares. 
Providing  that  the  matrix  C  is  not  singular,  these  equations  can  always  be  solved,  and 
the  formal  solution  xm  may  be  written 

xm  =  -C_1c0 

At  the  minimum,  the  merit  function  becomes 

<t>2  =  a0  •  a0  4-  c0  •  xm 

In  fact,  the  matrix  C  is  frequently  nearly  singular,  so  that  direct  matrix  inversion  is  an 
inappropriate  numeric  method  for  solving  the  normal  equations.  Before  exploring  a  more 
appropriate  method,  we  present  an  alternative  formulation  of  the  matrix  equations,  using 
homogeneous  matrices,  and  discuss  some  problems  associated  with  nonlinear  effects. 
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3.2.1  Homogeneous  Matrix  Formulation 


Homogeneous,  or  augmented  matrices  A *  and  x\  can  be  constructed  such  that 

f  =  a0  +  Ax  =  A'x'  =  (  A  |  a0  ) 

where  A'  has  N'  =  N  +  1  columns,  and  x'  has  N1  rows.  The  augmented  form  of  the 
covariance  array  C  is  given  by: 


;'  =  (- 


(c 

Co  \ 

[  4 

a0  '  a0  / 

Homogeneous  matrices  provide  a  compact  way  of  representing  and  storing  the  compo¬ 
nents  of  the  normal  equations, 

C'x'  =  0 


Since  C'  is  a  symmetric  matrix,  a  triangular  partition  can  be  used  to  represent  the 
full  matrix.  The  homogeneous  C  matrix  may  be  multiplied  by  a  nonzero  scalar  value 
without  changing  the  solution  of  the  system  of  equations.  For  example,  we  could  choose 
to  normalize  the  matrix  so  that  the  lower  right  corner  element  is  unity.  This  element 
is  zero  only  if  the  current  design  represents  an  ideal  solution,  in  which  case,  no  further 
processing  is  necessary. 


3.2.2  Nonlinear  Effects 

A  Taylor  series  expansion  of  a  scalar  function  of  several  variables  may  be  written  as 

<t>2  =  <t>l  +  go  ■  x  +  ^x*Hx 

where  H  is  the  Hessian  matrix  of  second  derivatives,  defined  as 


The  Hessian  matrix  may  be  expressed  as 


A  linear  defect  model  predicts  that 


4>2  =  a0  •  a0  +  2c0  •  x  +  x'Cx 

and  matches  the  first  and  second  terms  of  the  Taylor  series  but  it  leaves  out  the  second 
derivatives  in  the  Hessian  matrix. 


3.2.3  Damped  Least-Squares 

As  an  alternative  to  performing  lengthy  calculations  of  second  derivatives,  a  damping 
term  of  the  form  pi,  where  I  is  the  identity  matrix,  may  be  incorporated  to  eliminate 
singularities  and  restrict  the  least-squares  solution  to 

x  =  -  (C  +  pi)-1  c0 

In  the  limit  as  p  — >  0,  the  solution  approaches  the  least-squares  prediction.  An  exact 
solution  is  obtained  if  the  defects  are  linear  functions  of  the  variables.  Sufficiently  large 
values  of  p  guarantee  a  nonsingular  matrix  inversion  by  making  the  eigenvalues  non¬ 
negative.  As  p  becomes  large,  the  solution  approaches  a  small  step  opposed  to  the 
direction  of  the  local  gradient,  which  usually  results  in  a  smaller  merit  function  and  thus 
a  design  improvement.  The  method  of  damped  least-squares  has  been  used  extensively 
and  successfully  in  optical  design.  Damping  plays  an  important  role  in  any  nonlinear 
optimization  program,  but  we  choose  to  introduce  it  at  a  later  stage  in  the  process.  Let 
us  turn  our  attention  instead  to  the  technique  of  orthogonal  variables. 


3.3  Method  of  QU  Factorization 

QU  factorization  is  a  method  of  transforming  the  normal  equations  into  a  set  of  or¬ 
thonormal  equations.  We  seek  to  transform  the  matrix  A'  into  another  (M  x  N')  matrix 
Q  whose  columns  q,  are  orthonormal  over  the  defects, 

to 

<h  q*  =  n  QaQik  = 

t=l 
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The  constant  A;  represents  the  magnitude  of  q;,  It  is  zero  only  if  the  corresponding  q; 
vector  is  zero.  The  transformation  will  be  an  augmented  upper  triangular  matrix  U  such 
that 

QU  =  A' 

The  matrix  U  has  the  form 


1 

c*12 

£*13 

•  •  •  £*ijv 

£*io  ^ 

0 

1 

£*23 

•  *  ‘  £*2Af 

£*20 

0 

0 

1 

•  •  •  £*3N 

£*30 

0 

0 

0 

...  x 

£*ato 

Vo 

0 

0 

...  o 

i  / 

where 


q>  •  afc 

**  =  ~aT 


We  define  a  new  set  of  variables  w  such  that 

w  =  Ux' 

The  defect  vectors  can  now  be  written  as 

N 

f  =  A'x'  =  Qw  =  a0  +  <bwi 

j= i 

The  merit  function  itself  can  be  expressed  as  a  quadratic  sum, 

<f>2  =  +  E  +  «jo)2 

j= i 

where  the  minimum  value  of  <f> 2  is  be  obtained  by  inspection  as  w,  =  -aj0.  The  solution 
can  be  expressed  compactly  as 

w  =  -u0 

The  following  steps  are  repeated  as  often  as  necessary  to  obtain  a  final  solution: 

1.  Calculate  the  required  derivatives. 

2.  Carry  out  the  QU  factorization. 

3.  Determine  the  predicted  location  of  the  minimum  by  setting  w  =  — u0. 

4.  Solve  the  triangular  matrix  equation  w  =  Ux'  for  x. 

5.  Search  for  improvement  in  the  direction  predicted  by  the  model. 

Each  step  will  now  be  described  in  detail. 
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Calculation  of  Derivatives 

Derivatives  of  the  defect  functions  may  be  calculated  either  analytically  or  numerically. 
In  raytracing  parametric  surfaces,  the  partial  derivatives  are  calculated  analytically. 

These  derivatives  are  given  by 

Fu  =  S“  -  (Su • k)k 
Fv  =  Su  -  (Su  ■  k)k 


Cholesky  Decomposition 

The  Cholesky  decomposition  algorithm  is  one  way  of  carrying  out  the  factorization  of 
the  C  matrix  indicated  in  the  previous  section. 

The  general  column  qk  is  given  by 

k-l 

q*  =  ak  -  X] 

.=1 


where  for  j  <  k 


q,  •  qk  =  qy  afc  -  A*oc;k  =  0 


qj  ak 

a* = — r 


and  for  j  =  k 

The  Cholesky  algorithm  involves  the  following  a  recursive  process: 


Aj|  =  q*  •  q*  =  q*  • a* 


and 


a}k  ~  (a;  •  ak  [q»  '  aj] 


A*fc  -  C’fck  ~  Q>fcA> 

1=1 
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The  algorithm  is  a  simple  iterative  procedure.  For  the  first  column,  qi  =  at  with 
A l  =  Cn.  The  next  column  is  constructed  by  requiring  that 


Q2  —  a2  —  qi<*i2 


where 


<*12 

*2 


C*12 

A? 

C22  —  <*12  Aj 


Each  new  column  is  constructed  to  be  orthogonal  to  the  previous  columns.  This  is 
the  general  principle  of  Gram-Schmidt  orthogonalization  except  that  the  columns  are 
calculated  from  the  covariance  matrix  rather  than  the  original  A'  matrix.  If  the  sequence 
of  calculations  causes  a  negative  value  of  A2,  we  set  that  coefficient  to  zero  and  then 
continue  the  calculations.  If  the  merit  function  is  independent  of  a  particular  variable  Xj 
then  the  matrix  element  will  be  zero.  The  coefficients  0^  and  A2  must  also  be  zero. 
If  the  variable  is  linearly  dependent  on  earlier  variables,  then  the  algorithm  used  for  A2 
should  sum  to  zero  (or  within  a  specified  tolerance  of  zero). 


Obtaining  a  Solution 


The  merit  function  can  be  reduced  to  a  quadratic  form  through  the  following  steps: 

<j>2  =  f-f 

N  N  N 

=  a0  •  a0  +  2  £(a0  •  q>)«/,  +  £  £ (q,  •  qk)wjWk 

:=i  j—l t=i 

=  A2  +  £  A2(a20  +  2 aj0Wj  +  u >)) 

;=i 

=  A2  +  ^2  A2(t Vj  +  aj0)! 
i=i 


where 


A£  =  a0  •  a0  -  Y,  Ajajo 
i= 1 


The  minimum  value  of  <f> 2  can  be  obtained  by  inspection  as  w,  =  -aj0.  The  solution  can 
be  expressed  compactly  as 

w  =  —  u0 


27 


The  eigenvalues  {A2}  determine  the  relative  influence  of  the  transformed  variables  on  the 
magnitude  of  the  merit  function.  If  A 2  =  0,  then  w 3  has  no  effect  on  the  solution.  Zero, 
or  extremely  small  values  of  A2,  are  common  in  least-squares  equations.  Such  values 
results  from  variables  which  are  either  linearly  related  to  other  variables  or  which  have 
no  effect  on  the  merit  function.  This  situation  leads  to  an  ill-conditioned  matrix  inverse 
(indeterminate  values).  Using  QU  factorization,  we  can  recognize  variables  for  which 
A2  =  0,  and  exclude  them  from  further  processing.  Small  values  of  A2  are  associated  with 
vectors  oriented  along  valleys  of  the  merit  functions.  Larger  values  of  A2  are  associated 
with  directions  normal  to  the  walls  of  the  valleys.  If  A  is  zero,  we  have  an  unused  degree 
of  freedom  in  the  design.  There  is  a  continum  of  designs  that  satisfy  the  merit  function. 

Round-off  errors  can  lead  to  imaginary  values  for  some  A,  rather  than  small  positive 
values.  If  we  set  Ay  to  zero  at  this  stage  of  the  calculation,  we  can  avoid  some  of  the 
problems  of  ill-conditioning. 

Back  Substitution 
A  triangular  matrix  equation,  such  as 

w  =  Ux' 

is  easy  to  solve  for  x'  by  using  the  method  of  back  substitution.  First,  we  use  the 
definition  of  an  upper  triangular  matrix  to  write 

N 

Wi  =  Xi  +  J2  aHxi 
i=i+i 

Then  beginning  at  the  bottom  of  the  matrix,  we  can  generate  the  following  sequence: 

x/sr  =  tv/v 

Xff-1  =  Wjv-l  —  OtN-i'ffXN 

N 

xi  =  w>  ~  a'ixi 

j=i+i 

Search  in  One  Dimension 

Searches  in  one  dimension  may  be  made  along  any  single  variable  x},  transformed  variable 
Wj ,  or  solution  vector  u0.  The  general  case  is  to  define  a  search  vector  v  and  parameter 
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p  such  that 


x  =  pv 


where  p  =  0  is  the  current  point  and  p  =  1  is  the  predicted  solution  to 

<£{p)  -  0 

Our  technique  is  to  first  calculate  <£(1).  If  the  merit  function  is  smaller  than  the  original 
value  <£(0),  we  calculate  <£(1.5).  Otherwise,  we  find  <£(0.5).  At  this  stage,  we  have  enough 
information  to  construct  an  interpolating  parabola  <£(p)  =  50  +  bip  +  btp2.  We  locate 
either  the  minimum  value  or  the  closest  root  pm  and  then  find  <£(pm)  directly.  The  value 
of  p  associated  with  the  smallest  value  of  the  merit  function  is  used  to  generate  the 
starting  position  for  the  next  iteration  of  optimization. 
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Section  4 


Differential  Rays 


Differential  raytracing  is  a  simplified  method  of  finding  rays  that  are  close  to  some 
reference  ray.  If  the  reference  ray  is  the  central  ray  traveling  along  the  optic  axis,  then 
the  close  rays  are  called  paraxial  rays.  A  differential  raytrace  is  the  generalization  of  the 
idea  of  paraxial  rays.  Differential  rays  are  faster  to  calculate  than  exact  rays  and  can  be 
used  to  determine  the  focussing  properties  of  a  narrow  bundle  of  rays,  like  those  entering 
the  eye  from  an  external  object  point. 

Differential  rays  are  used  to  calculate  the  astigmatism  and  defocus  of  a  small  bundle 
of  rays.  Given  a  central  ray  from  p  traveling  along  the  unit  direction  k,  the  differential 
ray  is  given  by  6p  and  6k. 


4.1  Surface  Intersection 


The  differential  ray  intersects  a  surface  at  a  distance  q  +  6q  along  the  ray,  yielding  6 p' 
subject  to  the  condition  that  the  differential  ray  lies  in  the  surface  tangent  plane, 

Sp'  •  N  =  0. 

Taking  the  derivative  of  the  ray  translation  equation  gives 

6p'  =  6p  +  q  6k  +  5q  k 


Let 


p  =  6p  +  q  6k. 


Then 

and 


6q=  - 


p-N 

k-N 


6p'  =  p  4-  6q  k. 
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4.2  Ray  Refraction 


The  refracted  ray  is  given  by 
and  the  differential  ray  is 


n'k'  =  nk  +  *yN 


n'  <5k'  =  n  6k  +  7  <5N  +  67  N 
subject  to  the  condition  that  k'  is  a  unit  vector,  so  that 


k'  Sk'  =  0 


Let 

Then 

and 


k  =  n  £k  +  7  6N 


x  k'  k 

1  k'-N 


ri  $k'  =  k  +  67  N 


4.3  Differential  Trace  of  Parametric  Surface 


For  a  parametric  surface,  the  ray  intersection  differential  S S  is  given  by 

6S  =  S“  Su  +  Su  Sv 


Then 


Su  Su  4-  Sw  Sv  =  p  +  k  Sq 

represents  a  set  of  three  linear  equations  in  three  unknowns  (Su,  Sv,  Sq)  which  may  be 
solved  using  Cramer’s  rule: 


Su 

Sv 

Sq 


k  •  p  x  Su 

k  N 
k  •  Su  x  p 

k  •  N 
p  •  N 

~  k  •  N 
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The  differential  surface  normal  <*iN  is  also  required.  The  surface  normal  for  parametric 
surfaces  is  given  by 

N  =  S“  x  S" 

from  which  the  differential  normal  can  be  obtained  as 

6N  =  (<5S°)  x  S"  +  Su  x  (6SU) 

where 

SSU  =  Suu  Su  +  Suv  Sv 
SSV  =  Suu  Su  +  Svv  Sv. 
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Section  5 


Anisotropic  Refractive  Index  Ellipsoid 


This  section  explains  how  OPTRAN  computes  the  coefficients  of  an  ellipsoid  whose 
principal  axes  define  the  orthotropic  indices  of  refraction  at  a  particular  location  as  a 
function  of  temperature  and  stress. 


5.1  Orthotropic  Indices  of  Refraction 


The  equation  for  the  index  ellipsoid  of  an  unstressed  material  with  respect  to  the  principal 
axes  of  the  material  (x,  y,  z)  is  given  by 


x2  y2  z 2 
~r  +  ^  +  ~z  =  I- 

nl  nl  nl 


The  orthotropic  indices  of  refraction  (nx,  ny,  nz)  are  given  by 

dnx 


nx  =  n*  +  -rfAT 


n«  =  nu  + 


dT 

dnv 

dT 

dnr 


A  T 


n*  =  + 

where  (nz,  ny,  ht)  is  a  set  of  orthotropic  indices  of  refraction,  ( dnz/dT ,  dnv/dT,  dnz/dT) 
is  a  set  of  orthotropic  temperature  coefficients  along  the  same  axes,  and  AT  is  the 
temperature  change. 

A  geometric  transformation  is  required  if  the  axes  of  the  material  are  not  aligned 
with  the  global  axes. 
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5.2  Geometric  Transformation 


The  direction  cosines  of  the  material  axes  are  calculated  as  direction  cosines  with  respect 
to  a  set  of  global  coordinates  and  stored  as  column  vectors  in  the  matrix  X.  If  the  axes 
of  the  material  are  aligned  with  the  global  axes,  this  matrix  is  the  unit  matrix.  A  user 
provided  transformation  matrix  may  be  defined  for  materials  whose  axes  have  a  fixed 
orientation  with  respect  to  a  set  of  reference  axes,  which  may  be  the  global  axes  or  an 
orthogonal  set  of  axes  aligned  with  respect  to  the  parametric  parameters.  This  option 
was  intended  to  accommodate  materials  that  are  deformed  to  follow  an  irregular  shape. 

The  choices  are  distinguished  by  setting  the  following  flags  for  each  material: 


IUVWAX(IMATL)  =  0  use  the  global  axes  as  reference  axes 

■  1  use  the  derivatives  of  the  parametric  parameters 
with  respect  to  the  global  axes  to  define  a  set 
of  reference  axes 


IOPMTR(IMATL) 


0  do  not  use  the  OPMATR  matrix  to  define  the  material 
axes  with  respect  to  the  reference  axes 
1  use  the  OPMATR  matrix 


Let  D0  contain  the  direction  cosines  of  the  material  axes  with  respect  to  the  reference 
axes  stored  as  column  vectors, 


Do  =  [  J  m  n  ] 

where  I,  m,  n,  are  unit  direction  cosine  vectors  which  define  the  direction  of  the  three 
material  axes  with  respect  to  the  reference  axes. 

If  the  global  axes  are  used  as  the  reference  axes  X  is  simply  set  equal  to  input  material 
axis  direction  cosine  matrix  OPMATR  for  material  IMATL. 


X  =  Do 


Otherwise  a  set  of  orthogonal  reference  axes  is  computed  as  a  function  of  the  derivates 
of  the  x,  y,  z  coordinates  with  respect  to  the  u,  v,  w  geometric  hyperpatch  parametric 
parameters.  Then 

X  =  DoDuuu, 
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where 


Duuu,  =  [vt  V2  V3  ] 


The  vectors  defined  by  partial  derivative  triples  Vu  =  ( Xu ,  Yu,  Zu),  Vv  =  ( Xv ,  Yv, 
Zv),  Vw  —  {Xw,  Yw ,  Zw)  are  not  orthogonal  except  in  special  cases  where  the  geometry 
has  curvilinear  orthogonality.  An  orthogonal  set  of  direction  cosines  is  computed  by 
taking  cross  products  of  partial  derivative  triples  and  normalizing  the  results  to  unit 
vectors.  In  general  only  the  direction  of  one  of  the  vectors  remains  unaltered  and  can  be 
used  as  reference.  By  an  arbitrary  decision,  the  following  convention  has  been  adopted 
for  defining  the  a  reference  axis  system  (V|,  V2,  V3)  based  on  the  parametric  partial 
derivatives  is  computed  as  follows: 

If  IVOPT  =  1  (  Parameter  u  held  constant  ) 

v2  =  Vj\Vv\ 

Vi  =  (VvxVw)/\V„xVw\ 

V3  =  (Vi  x  Vt)/\Vi  x  Vt\ 

If  IVOPT  =  2  (  Parameter  v  held  constant  ) 

Vi  =  VJ\VU\ 

v2  =  (VWXVU)/ \VmxVu\ 

V's  =  (V,  x  Vjj/IV,  x  V2\ 

If  IVOPT  =  3  (  Parameter  w  held  constant  ) 

Vi  =  VJ\VU\ 

v>  =  (V.  x  V.)i\V,  x  V,| 

V,  =  (V,  x  Vt)l\Vt  x  V,| 

where  the  program  variable  IVOPT  is  set  as  a  function  of  the  incident  face. 


5.3  Stress  Birefringence 

When  a  stress  is  applied  to  a  material,  the  index  ellipsoid  is  modified,  and  the  changes  in 
the  components  of  the  dielectric  tensor  are  linearly  related  to  the  six  stress  components. 
[3|(p.  703-704] 
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On  applying  stress,  the  index  ellipsoid  is  changed  into  another  whose  equation  is 
+  avvy 2  +  a«22  4-  2  axvxy  +  2  avzyz  4-  2axzxz  =  1. 

The  coefficients  of  the  ellipse  are  defined  by 

{A}  =  {A0}  +  \q]{a} 

where 


The  first  term  contains  the  orthotropic  indices  and  associated  thermal  coefficients.  The 
second  term  contains  stress-related  terms  where  {a}  is  the  stress  tensor  and  [g)  is  a  set 
of  36  stress-optical  coefficients. 

Although  all  the  stress-optic  coefficients  may  have  unique  nonzero  values,  symmetry 
considerations  for  crystalline  or  isotropic  materials  prescribe  relationships  among  the 
coefficients  and  reduce  the  number  of  independent  values  that  must  be  specified.  For 
cubic  crystals,  the  three  principal  axes  (x,  y,  z)  are  equivalent,  and  consequently  the 
following  relations  hold  among  the  stress-optics  coefficients: 

9n  =  <722  =  <733 

<712  =  <721  =  <723  =  *732  =  <7l3  =  <731 

<744  =  955  =  966 

with  all  the  remaining  coefficients  being  zero. 

For  isotropic  materials,  the  above  relations  must  remain  unaltered  for  any  change  of 
axes.  This  is  only  possible  if  the  stress-optical  coefficients  satisfy  the  additional  relation 

2 9<4  =  9u  -  9i2- 

If  the  global  coordinate  system  is  different  from  that  of  the  dielectric  material,  the 
stresses  must  first  be  converted  to  the  material  coordinates.  Then  the  resulting  dielectric 
tensor  must  be  transformed  back  into  the  global  coordinate  system.  If  {<?}'  denotes 
material  coordinates  and  {cr}  global  coordinates,  then 

W'  =  NW 
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where  (sj  is  the  coordinate  transformation  matrix.  In  material  coordinates, 

{A}'  =  {Aoy  +  [g)'{ay 

which  may  be  written  in  reference  coordinates  as 

{A}  =  {s\~l{Aoy  +  {sr\q\'[s\{a}. 

The  matrix  [<s]_1[g]'[s]  is  not  a  function  of  stress  and  in  many  cases  is  not  a  function  of 
location,  and  therefore  can  be  precomputed. 


5.4  Tensor  Rotation  Transformation 

Given  a  point  expressed  in  (z,  y,  z)  global  coordinates,  we  find  the  corresponding  repre¬ 
sentation  in  (£,  m,  n)  coordinates.  Let  1,  m,  n  be  the  unit  direction  of  the  material  axes 
in  (z,  y,  z)  global  coordinates.  Then 


(  £  m  n)  =  (z  y  2 ) 


L  m, 


n. 


£*  m, 


my  riy 


n2 


The  reverse  transformation  is 


(*  y  *] 

1  = 

£  m  n  ) 

'  £* 
mx 

£y 

mv 

l,  ' 
mz 

.  nz 

nv 

n*  . 

The  coordinate  transformation  matrix  [sj  is  given  by 


a 

i2 

2tXiy 

2  tyix 

2  £,£. 

ml 

m] 

2  mxmv 

2  mvmx 

2  mxmx 

ni 

nl 

< 

2nxnv 

2  nvnz 

2  Tlx7lx 

f-zmx 

tvmv 

izmt 

(£xmy  +  mx£v) 

(tvmz  +  mz£y) 

(£lml  4-  £xttix) 

mtnz 

mvnv 

mxnx 

( mxnv  +  nxmv) 

(mvnz  4-  nxmv) 

(mxnx  4-  m,nx) 

txTlj 

[tznv  4-  nzlv) 

(£„n,  4-  nzly) 

(£r^r  4" 

37 


5.5  Principal  Axes  of  Projected  Ellipse 

Given  an  ellipsoid  A  in  the  (x,  y,  z)  coordinate  system  and  an  axis  of  projection  given 
by  the  unit  vector  (u,  v,  w).  Then 


ux  4  vy  4  wz  —  0 


is  the  equation  of  a  plane  through  the  origin  and  perpendicular  to  the  projection  axis. 
The  intersection  of  this  plane  with  the  ellipsoid  is  an  ellipse.  We  want  to  find  the  major 
and  minor  axes  of  this  ellipse. 

We  find  the  lengths  of  the  axes  by  finding  the  extrema  of  r2  =  x2  4  y2  4  z2  subject 
to  the  constraints  that  the  extrema  lies  on  the  plane  and  the  ellipsoid. 

By  using  Lagrange’s  method  of  undetermined  multipliers  (Ai  and  A2),  we  define  the 
following: 

h{x,y,z)  =  x2  4  y2  4  z2  - (-  A ihx{x,y,z)  4  A2h2(x,y,z) 

where 


hi(x ,  y,  z)  ~  axxx 2  4  avvy2  4  aztz 2  +  2 axyxy  4  2 ayzyz  4  2 axxxz  -  1 
/i2(x,y,  z)  =  ux  +  vy  +  wz. 

Then  setting  the  partial  derivatives  of  h(x,  y,z)  =0  gives 
dh 


dx 

dh 

dy 

dh 

dz 


=  2x  +  2Ax  [a  xxX  -h  Uxyy  "t"  UyzZ J  -f-  A2u  —  0 
=  2y  +  2Aj  [ovyy  4-  axyx  +  ayxz ]  +  A2v  =  0 
=  2z  4  2Ai  \atzz  4  axxz  4  ayzy\  4  A 2w  =  0. 


Then  we  may  obtain 


and 


where 


dh  dh  dh  , 

x—  +  y-z-  +  z—  =  2 r2  4  2Ai  =  0 
dx  dy  dz 


dh  dh  dh  t  .  . 

u- — I-  v—  4  —  =  2Ai  [uax  4  vay  4  vuaz j  4  A2  =  0 


dx 


dy 


ax  —  ctxxx  4  4  uxzz 

uy  —  uTyX  4  avyy  4  uyzz 

uz  —  uxzx  4  Qyzy  4  uzzz> 
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Then 


=  2x  -  2 r2ax  4-  2 r2u  j uax  +  vay  4-  tua2]  =  0 

=  2 y  -  2 r2av  4-  2 r2v  [uat  4-  vay  4-  tya2]  =  0 

=  2z  —  2 r2az  +  2r2tt>  \uax  4-  vay  4-  u>az j  =  0. 
Factoring  gives  the  following  system  of  homogeneous  equations. 

x[r~2  -  axx  4-  uau\  +  y[-axy  4-  uav ]  4-  z[-axx  4-  uaw\  = 
x\-axy  +  t>au]  4-  y[r~2  —  ayy  +  va„]  +  z[-ayz  4-  uaw]  = 
x[-dxz  +  u>au]  -h  y[-ayz  4-  wav\  4-  z\r~ 2  -  azz  +  waw ]  = 

where 

a„  =  +  dxyv  +  ax*w 

tty  —  dxyu  +  ayyv  4*  dyzxv 
aw  =  axzu  4-  ayz v  4-  azzw. 


dh 

dx 

dh 

dy 

dh 

dz 


0 

0 

0 


The  determinant  of  coefficients  of  this  set  of  three  homogeneous  equations  must  vanish 
for  there  to  be  a  nontrivial  solution  for  x,  y,  and  2.  Using  this  requirement  gives 


“2  -  dxx  4-  uau] 

[  dxy  +  ua„] 

\-axz  +  udw\ 

[  dxy  4“  vUu] 

[r~2  -  dyy  +  va„j 

[-a„z  4-  va,u] 

\-a-xz  4-  umu] 

[-Oy*  +  waA 

[r~2  -  azt  4-  tvd , 

Expanding  the  determinant  gives 

ar4  +  6r2  4- 1  =  0 


where 


a  —  (flj  —  o-i)  —  (03  —  d4) 


fll  —  dxxdyy  4"  dyydzz  4"  dxxdzz 

=  a2„  +  aj,  + 

03  —  dxx{vdv  4~  d) j  4"  4“  4"  0Lzz\xto>^i  4~ 

d4  =  av2(vaw  4-  wav)  4-  axy{vau  4-  ua„)  4-  axz(uaw  4-  wau), 

and 

b  —  Ufly  4"  Vdfj  4“  tVdyj  dXX  flyy 

This  equation  is  a  simple  quadratic  in  r2.  The  solutions  are  the  major  and  minor  axes. 
The  system  of  homogeneous  equations  can  then  be  solved  for  the  direction  of  the  respec¬ 
tive  principal  axes. 


Section  6 


Optical  Waves  in  Anisotropic  Materials 


Transparent  materials  may  be  optically  anisotropic  due  to  intrinsic  crystal  properties 
or  as  a  result  of  mechanical  stresses.  In  either  case,  the  index  of  refraction  becomes  a 
function  of  direction  within  the  material  and  can  be  described  by  the  equation  of  an 
ellipsoid.  An  anisotropic  material  permits  two  monochromatic  plane  waves  with  two 
different  polarizations  and  two  different  velocities  to  propagate  in  any  given  direction. 
The  directions  of  the  two  displacement  field  vectors  D  corresponding  to  a  given  direction 
of  propagation  k  are  perpendicular  to  each  other.  The  phenomonon  of  wave  propagation 
in  an  anisotropic  material  is  called  birefringence  or  double  refraction. 


6.1  Light  Propagation  in  Anisotropic  Materials 


As  discussed  in  the  previous  section,  the  equation  for  the  index  ellipsoid  with  respect  to 
the  principal  axes  of  the  material  (z,  y,  z)  is  given  by 


=  1. 


After  rotating  the  coordinate  system,  this  ellipsoid  is  changed  into  another  whose  equa¬ 
tion  is 

oZIz  +  -i-  “f*  2fljyZy  -t-  2oyZyz  "f-  2ozzzz  1. 


In  a  birefringent  medium  two  orthogonal  polarized  waves  may  be  propagated  in  a 
given  direction.  Let  a  set  of  principal  axes  (z,  y,  z )  be  defined,  with  the  direction  of 
propagation  in  the  z-direction,  the  allowed  polarization  directions  along  the  x  and  y  axes, 
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Figure  6.1:  Construction  of  the  D  vectors  belonging  to  a  wave  normal  k  [lj. 


and  the  associated  indices  of  refraction  be  given  by  nx  and  nv.  The  Jones  matrix  for  the 
propagating  wave  is  then 


J'  = 


e-2  rjnxd/X 


0  e~2irjn,d/\  j  *» 

where  d  is  the  distance  propagated.  This  relationship  may  be  factored  by  defining 


n  = 


nz  +  n„ 


n,  —  ny 
0n  —  - r - • 


The  optical  path  S  along  the  ray  is  given  by 

6  =  nd 


and  the  Jones  matrix  is  now  given  by 


J'  = 


/  e-2x}6nd/\ 

\  0 


0 

e2)r  jSnd/)> 


) 


J. 


which  may  be  expressed  as  the  Mueller  matrix  of  a  linear  retarder. 


As  shown  in  Figure  6.1,  the  intersection  of  the  dielectric  ellipsoid  with  a  plane  through 
the  origin  and  perpendicular  to  the  direction  of  wave  propagation  k  defines  an  ellipse 
whose  major  and  minor  axes  determine  the  direction  of  the  displacement  field  vectors  d 
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and  the  associated  indices  of  refraction  for  the  two  orthogonal  polarized  waves  propagat¬ 
ing  in  the  direction  k.  The  vector  d  is  a  unit  vector  in  the  direction  of  the  displacement 

D. 

Wavefront  propagation  is  governed  by  the  triplet  of  unit  vectors  (d,  h,  k)  such  that 

k  =  d  x  h 

Rays  in  a  birfringent  medium  do  not  propagate  in  the  same  direction  as  the  wavefront. 
Ray  propagation  is  determined  by  the  Poynting  vector  S  given  by 

S  =  E  x  H 

The  direction  of  ray  propagation  and  polarization  is  represented  by  the  triplet  of  unit 
vectors  (e,  h,  s)  such  that 

s  =  exh 

The  electric  field  E  is  related  to  the  displacement  field  D  by  the  dielectric  tensor. 

D  =  eE 

The  index  ellipsoid  is  an  equivalent  representation  of  the  dielectric  tensor  that  relates 
the  electric  and  displacement  fields  as  follows: 

/  Ez  ^  dzz  aZy  azx  l  Dz 

I  Ey  —  QZy  dyy  dy  X  I  Dy  (6.l) 

V  ^x  )  .  .  v  D z  ) 

The  indices  of  refraction  may  be  graphed  on  a  polar  plot  as  a  function  of  the  direction 
of  wave  propagation.  This  defines  a  two-sheeted  surface  called  the  normal  surface.  Where 
the  sheets  intersect,  the  two  orthogonal  polarized  waves  have  the  same  index  of  refraction. 
These  directions  are  called  the  optic  axes  of  the  medium. 

A  material  is  classified  as  isotropic  if  the  three  principal  indices  of  refraction  are 
equal,  as  uniaxial  if  two  of  the  three  principal  indices  are  equal,  and  as  biaxial  if  none 
of  the  principal  indices  are  the  same.  Table  6.1  lists  typical  refractive  indices  of  some 
crystals. 

The  wave  normal  surface  for  an  isotropic  material  is  a  sphere,  since  the  index  of 
refraction  does  not  vary  with  direction.  For  a  uniaxial  material,  the  normal  surface 
consists  of  a  sphere  and  an  ellipsoid  of  revolution.  If  nz  —  nv  then  these  two  sheets  touch 
at  two  points  on  the  z-axis  because  both  wavefronts  propagate  along  the  z-axis  with  the 
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Table  6.1:  Typical  refractive  indices  of  some  crystals  [3]. 


Crystal 

Refractive  Indices 

Isotropic 

Fluorite 

1.392 

Sodium  chloride,  NaCl 

1.544 

Diamond,  C 

2.417 

CdTe 

2.69 

GaAs 

3.40 

", 

Uniaxial 

Ice.  HjO 

1.309 

1.310 

(Positive) 

Quartz,  SiOj 

1  544 

1  553 

Beryllium  oxide,  BeO 

1.717 

1.732 

Zircon,  ZrSiO, 

1.923 

1  968 

ZnS 

2.354 

2.358 

Rutile,  TiOj 

2.616 

2.903 

", 

Uniaxial 

ADP,  (NH4)H,P04 

1.522 

1.478 

(Negative) 

Beryl,  BejAlj(SiOj)* 

1.598 

1.590 

KDP,  KH2P04 

1  507 

1.467 

Sodium  nitrate,  NaNO, 

1.587 

1.366 

Calcite,  CaCO, 

1.658 

1.486 

Tourmaline 

1.638 

1.618 

Sapphire.  AljO, 

1.768 

1.760 

Lithium  niobate  LiNbO, 

2.300 

2.208 

Barium  titanate,  BaTiO, 

2.416 

2.364 

Proustite,  AgjAsS, 

3.019 

2.739 

", 

", 

n. 

Biaxial 

Gypsum 

1.520 

1.523 

1  530 

Feldspar 

1.522 

1.526 

1.530 

Mica 

1.552 

1.582 

1.588 

Topaz 

1.619 

1.620 

1.627 

Sodium  nitrite 

1.344 

1.41 1 

1.651 

YAlOj 

1.923 

1.938 

1  947 

SbSI 

2.7 

3.2 

3.8 
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<*> 


(b) 


<c) 


Figure  6.2:  Intersection  of  the  normal  surface  with  xz  plane  for  (a)  biaxial  crystals,  (b) 
positive  uniaxial  crystals,  and  (c)  negative  uniaxial  crystals  [3 j. 


Figure  6.3:  Orientation  of  rays  and  waves  in  a  uniaxial  crystal  [2]. 
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same  index  of  refraction.  This  direction  is  called  the  optic  axis  of  the  material.  A  biaxial 
material  has  two  separate  optic  axes  coplanar  with  longest  and  shortest  principal  axes. 
Figure  6.2  shows  cross-sections  of  the  normal  surface  for  biaxial  and  uniaxial  crystals. 

Figure  6.3  shows  Huygens  wavelets  for  two  modes  of  propagation  in  an  anisotropic 
crystal.  Huygens  wavelets  for  which  the  index  of  refraction  is  independent  of  direction  are 
spherical  wavefronts,  and  the  associated  rays  are  called  ordinary  waves.  The  direction 
of  wavefront  propagation  and  the  direction  of  energy  transfer  (ray  direction)  are  the 
same.  Wavelets  for  which  the  index  of  refraction  varies  with  direction  are  ellipsoidal 
wavefronts,  and  the  associated  waves  are  called  extraordinary  waves.  The  direction  of 
ray  propagation  is  different  from  that  of  wavefront  propagation.  The  envelope  of  all 
wavelets  in  both  cases  is  a  plane  wave  propagating  horizontally  to  the  right. 

Given  the  direction  of  wave  propagation  k  and  the  refractive  index  ellipsoid  A,  we 
can  calculate  dlt  n1?  d2,  and  n2  for  the  principal  polarization  states,  as  shown  in  Section 
5.  If  ni  =  n2  there  is  no  unique  choice  of  principal  directions,  and  we  may  arbitarily 
select  the  principal  polarization  states. 

For  each  polarization  mode,  we  calculate  h  =  k  x  d  and  the  electric  field.  Then  using 
the  directions  of  the  electric  and  magnetic  fields,  we  can  calculate  the  ray  direction  s 
from 

s  =  e  x  h 

The  ray  index  of  refraction  nr  is  obtained  from  the  wave  index  of  refraction  by 

nr  =  nk  •  s 

For  a  wavefront  propagating  normal  to  a  plane-parallel  slab  of  material  of  thickness  d, 
as  shown  in  Figure  6.4,  the  optical  path  6  is  given  by 

6  =  nd 

If  the  ray  direction  makes  an  angle  a  with  the  wave  direction,  the  distance  traveled  by 
the  ray  in  the  material  is  d/  cos  a  and  the  index  is  ncosa,  so  that  the  calculated  optical 
path  is  the  same  whether  calculated  for  the  wavefront  or  the  ray. 

Figure  6.5  is  an  illustration  of  the  propagation  of  light  through  a  calcite  crystal.  If 
we  send  a  narrow  beam  of  normal  (unpolarized)  light  into  a  calcite  crystal,  it  will  split 
into  two  beams.  Rotating  the  crystal  causes  one  of  the  rays  to  remain  stationary  and 
the  other  to  move  in  a  circle  about  it,  following  the  motion  of  the  crystal.  The  fixed 


Figure  6.5:  A  light  beam  with  two  orthogonal  field  components  traversing  a  calcite 
principal  section  [2j. 
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Figure  6.6:  Double  refraction  at  the  boundary  of  an  anisotropic  material. 

ray,  which  follows  the  usual  rules  of  refraction,  is  the  ordinary  ray.  The  moving  ray,  not 
following  the  usual  rules,  is  the  extraordinary  ray. 

In  the  discussion  so  far,  the  incident  wavefront  has  been  normal  to  the  crystal  surface, 
so  that  the  direction  of  refraction  was  also  normal  to  the  crystal  surface.  In  general, 
however,  rays  will  be  incident  at  an  oblique  angle  to  a  dielectric  boundary  surface. 


6.2  Double  Refraction 


Consider  an  unpolarized  plane  wave  incident  on  the  surface  of  an  anisotropic  material. 
The  refracted  wave,  in  general,  is  a  mixture  of  two  propagation  modes.  The  conditions 
for  refraction,  established  by  the  application  of  Fermat’s  Principle,  required  that  the 
refracted  wavefront  lie  in  the  plane  of  incidence  established  by  the  incident  propagation 
vector  and  the  surface  normal  and  that  Snell’s  Law  is  satisfied, 

n  sin  0  =  rii  sin  9\  —  n?  sin  0? 

For  an  extraordinary  ray,  the  index  of  refraction  is  a  function  of  the  angle  of  incidence, 
leading  to  the  following  transcendental  equation 

nsin0  =  nr(0r)sin0r 

which  must  be  solved  for  0r.  The  solution  can  be  obtained  graphically  as  shown  in 
Figure  6,6.  The  wavefront  optical  cosines  are  vectors  whose  magnitude  is  the  index  of 
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Figure  6.7:  Ordinary  reflection  at  the  boundary  of  an  anisotropic  material. 

refraction.  The  wavefront  normal  surface  traces  the  index  of  refraction  as  a  function 
of  angle  of  incidence  for  the  two  possible  propagation  modes.  The  projection  of  the 
wavefront  propagation  vector  onto  the  boundary  is  equal  to  nsin0.  This  projection  must 
be  the  same  for  the  incident  and  refracted  rays  to  satisfy  Snell’s  Law.  The  refracted  waves 
may  be  found  by  extending  a  vertical  line  from  the  incident  wavefront  propagtion  vector 
through  the  normal  surface.  The  intersection  points  with  the  normal  surface  represent 
the  two  refracted  waves. 

If  a  polarized  wave  propagating  within  an  anisotropic  material  is  incident  at  an  exit 
surface,  part  of  the  wave  will  be  reflected  at  the  surface.  For  an  ordinary  wave,  shown 
in  Figure  6.7,  the  angle  of  reflection  will  be  equal  to  the  angle  of  incidence.  For  an 
extraordinary  wave,  shown  in  Figure  6.8,  the  angle  of  reflection  is  not  equal  to  the  angle 
of  incidence  because  the  index  of  refraction  is  different  for  the  reflected  wave. 

In  examining  Figures  6.7  and  6.8,  we  observe  that  there  are  two  possible  solutions 
for  the  reflected  wave.  Only  one  solution  is  reported  for  each  situation.  The  polarization 
modes  for  this  geometry  turn  out  to  be  uncoupled.  The  ordinary  wave  reflects  only  as 
an  ordinary  wave,  and  the  extraordinary  incident  wave  reflects  only  as  an  extraordinary 
wave.  This  situation  is  shown  in  Figure  6.9.  The  short  rays  denote  the  Poynting  vector 
and  the  longer  waves  the  wave  propagation  vector.  The  p-polarization  case  corresponds 
to  the  extraordinary  wave.  The  associated  wave  normal  surfaces  are  shown  for  each 
polarization  mode. 

In  general,  at  the  boundary  between  two  different  anisotropic  materials,  there  are 
two  possible  reflected  waves  and  two  refracted  waves  for  a  polarized  incident  wave.  The 
energy  carried  by  the  incident  wave  is  divided  among  the  four  possible  output  waves. 
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Figure  6.8:  Extraordinary  reflection  at  the  boundary  of  an  anisotropic  material. 

The  distribution  of  energy  is  determined  by  using  the  boundary  conditions  imposed  by 
Maxwell’s  equations  on  the  transverse  components  of  the  electric  and  magnetic  fields. 
This  distribution  will  be  discussed  in  the  next  section. 


6.3  Reflection  and  Refraction  Coefficients 


Let  s<,  be  the  unit  Poynting  vector  of  the  incident  ray,  e„  the  electric  field  vector,  and  h* 
the  magnetic  field  vector,  such  that 


So  —  ®o  ^  ho 

Let  z  be  the  normal  to  the  surface,  with  orthogonal  vectors  x  and  y  lying  in  the  plane 
of  the  surface.  There  will  be  two  reflected  beams,  sa  and  s,.  and  two  refracted  beams  sk 
and  s d.  At  the  boundary,  the  transverse  components  of  electric  and  magnetic  fields  are 
continuous. 

These  boundary  conditions  are: 

E0eo  ■  x  4-  Eaea  •  x  4-  Ecee  ■  x  =  Ebeh  •  x  4-  Eded  •  x 
E0e0  •  y  4-  Eaea  ■  y  +  Eeee  •  y  =  Ebeh  •  y  4-  Eded  •  y 
n0E0h0  ■  x  +  naEaha  ■  x  +  neEe he  •  x  =  n»£Vh»  •  x  4-  ndEd hd  •  x 
n0E0h0  •  y  +  naEaha  ■  y  +  ncEc hc  •  y  =  njE’jh*  •  y  +  ndEdhd  •  y 
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P-polarizatlon 


S -polarization 


Figure  6.9:  Refraction  and  reflection  at  boundary  between  two  anisotropic  materials. 
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-ea  x  e6  x  -ec  x  ed  x 
-ea  •  y  eb  ■  y  -ec  •  y  •  y 

-nahax  n6hi,-x  -nehe-x  ndhd  ■  x 
.~naha  y  nbhb  ■  y  -nehcy  ndhd-y 

Let  x  be  perpendicular  to  the  plane  of  incidence,  then 

X  =  S0  X  z 
y  =  z  x  x 


Ea  ' 

e0  x 

Eb 

G<>  •  y 

Ec 

n0 h0  •  x 

.  Ed  . 

.  ftoho • y 

If  the  Poynting  vector  is  parallel  to  the  surface  normal  (normal  incidence),  no  plane  of 
incidence  is  defined.  The  choice  of  x  is  then  arbitrary.  We  typically  let  x  =  e0. 

The  following  amplitude  reflection/transmission  coefficients  may  be  jtained 


where  m  is  one  of  the  four  reflecting  or  refracting  rays  (a,  b,  c,  d).  Energy  reflec¬ 
tion/transmission  coefficients  are  obtained  from 


Tm  = 


ftmSm 

n0 • 


•  z 
z 


2 


6.3.1  Reflection  and  Refraction  in  Isotropic  Media 

The  equations  for  reflection  and  refraction  at  the  boundary  between  two  isotropic  ma¬ 
terials  define  the  Fresnel  coefficients.  In  this  section  we  show  how  they  may  be  derived 
from  the  general  expressions  given  above. 


S-Polarization 

For  an  incident  ray  polarized  perpendicular  to  the  plane  of  incidence  (s-polarization), 
e„  =  x  and 


’  -1 

1 

0 

0  ' 

r  Ea 

1 

0 

0 

cos  8 

cos  8' 

Eb 

0 

0 

0 

na 

~n, 

Ec 

0 

n0  cos  8 

n,  cos  $' 

0 

0  . 

Ed  J 

n0cos8  _ 
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from  which  we  obtain  Ec  =  Ej  =  0  and 


-1  1 

'  Ea  ' 

1 

n0  cos  6  n,  cos  9' 

.  E»  . 

n0 cos  9 

These  equations  can  be  solved  for  the  reflection  and  transmission  amplitude  coefficients 

2 n0  cos  9 

n0  cos  9  +  n,  cos  9' 
n0  cos  9  -  n,  cos  9' 

'  n0  cos  9  -f  n,  cos  9' 


P-Polari*ation 


For  an  incident  ray  polarized  parallel  to  the  plane  of  incidence  (p-polarization),  e0  =  s0  xx 
and 


’  -1 

1 

0 

0  ' 

‘  Ea  ‘ 

'  0 

0 

0 

cos# 

cos#' 

Eh 

cos# 

0 

0 

n0 

-»• 

Ee 

-n„ 

nocos0 

n3  cos  0' 

0 

0 

Ed  J 

0 

from  which  we  obtain  Ea  =  Eb  —  0  and 


cos# 

n0 


cos  9' 
~n. 


'  Ec  ' 

cos# 

.  E*  . 

.  ~n°  . 

These  equations  can  be  solved  for  the  reflection  and  transmission  amplitude  coefficients 

Ed  2  n0  cos  9 

p  E0  n,  cos  9  +  n0  cos  9' 

Ee  n3  cos  9  —  n0  cos  9' 

p  E0  n,  cos  9  +  n„  cos  9' 


Normal  incidence 

For  isotropic  media  at  normal  incidence,  there  is  no  plane  of  incidence  defined.  We  can 
then  choose  x  =  e„.  Then  the  equations  simplify  to 

E0  +  E a  =  Eh 
naELl  -  naEa  =  nfEh 
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from  which  we  can  find  the  transmission  and  reflection  amplitude  coefficients 

^  __  Eb_  _  2  n„ 

E0  n0  +  n, 

^  _  Ea  Tls  TIq 

E0  n,  +  n0 
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